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Abstract 

QCD with three degenerate quark flavours at zero baryon density exhibits a 
first order thermal phase transition for small quark masses, which changes to a 
smooth crossover for some critical quark mass TOq, i.e. the chiral critical point. 
It is generally believed that as an (even) function of quark chemical potential, 
the critical point m c (/i) moves to larger quark masses, constituting the critical 
endpoint of a first order phase transition in theories with m > ra[j. To test this, 
we consider a Taylor expansion of m c (u) around /i = and determine the first two 
coefficients from lattice simulations with staggered fermions on N t = 4 lattices. 
We employ two different techniques: a) calculating the coefficients directly from 
a fj, = ensemble using a novel finite difference method, and b) fitting them 
to simulation data obtained for imaginary chemical potentials. The a 2 and fi 4 
coefficients are found to be negative by both methods, with consistent absolute 
values. Combining both methods gives evidence that also the u 6 coefficient is 
negative. Hence, on coarse Nt = 4 lattices a three-flavour theory with m > ttiq 
does not possess a chiral critical endpoint for quark chemical potentials u<T. 
Simulations on finer lattices are required for reliable continuum physics. Possible 
implications for the QCD phase diagram are discussed. 



1 Introduction 

The understanding of the QCD phase diagram as a function of temperature and baryon 
density plays an important role in an increasing number of nuclear and particle physics 
programmes, with possible applications to cosmology and astro-particle physics. Of 
particular interest in the context of heavy ion collision experiments is the location 
and nature of the quark-hadron transition. The generally expected picture - with 
an analytical crossover at /z = turning into a first order transition beyond some 
critical chemical potential /ig - is based on combining lattice results for /i = with 
model calculations at /i ^ and connecting various limiting cases by universality and 
continuity arguments pQ. In particular it is generally believed that, in the enlarged 
parameter space {m u d) m s , T, //}, the critical endpoint of QCD with physical quark 
masses is continuously connected to the chiral critical line at \Ie — 0, which delimits 
the quark mass region with a first order chiral phase transition at zero density. 

The schematic situation at fi = is depicted in Fig. [T] (left). In the limits of zero 
and infinite quark masses (lower left and upper right corners), order parameters cor- 
responding to the breaking of a global symmetry can be defined, and one numerically 
finds first-order phase transitions at small and large quark masses at some finite tem- 
peratures T c (m). On the other hand, one observes an analytic crossover at intermediate 
quark masses, with second order boundary lines separating these regions. Both lines 
have been shown to belong to the Z(2) universality class of the 3d Ising model [SHU [5]. 
The chiral critical line has been mapped out on N t = 4 lattices, and the physical point 
is confirmed to be on the crossover side of this line [HI E] ■ The continuity arguments in 
[1] suggest that this line should shift with /i towards larger quark masses, thus increas- 
ing the region of a chiral first order transition. At some value \xe it passes through 
the physical point, whence it represents the endpoint of the phase diagram for physical 



This scenario has recently been called into question by numerical simulations at 
imaginary chemical potential /i = /i; 6 1 [6J. There it was found that, as a function 
of fii, the chiral critical line moves towards larger quark masses. In the degenerate 
Nf = 3 theory and in the neighbourhood of the critical point m c ud = m c s = mg, the 
dependence of the critical temperature and quark mass on the chemical potential can 
be described by a Taylor expansion, with only even powers of /i due to CP-symmetry, 



QCD, Fig. HI (middle). 
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Figure 1: Left: Schematic phase transition behaviour of Nf = 2 + 1 QCD for different 
choices of quark masses (m u ^,m s ) at fi = 0. Middle/Right: Critical surface swept 
by the chiral critical line as fi is turned on. Depending on the curvature, a QCD 
chiral critical point is present or absent [2J. For heavy quarks the curvature has been 
determined [5] and the first order region shrinks with \x. 

By performing simulations at several values of /i, and quark masses, the leading coef- 
ficients were determined, and in particular c\ was found to be negative. This implies 
a scenario as in Fig. [1] (right), where the chiral critical line does not cross the physical 
point for small (real) chemical potentials. If this behaviour continues to larger values of 
/j,, there would be no non-analytic chiral phase transition or end point at all in physical 
QCD. 

Clearly, it is now important to assess the reliability of this result. The main sys- 
tematic uncertainties associated with the calculation in [6] are the truncation of the 
Taylor expansion Eq. (|2D and the coarse lattice with Nf = 4, a ~ 0.3 fm. In this letter 
we perform a systematic investigation of the former. Firstly, instead of extracting the 
leading term in Eq. ([2j) from fits to a truncated Taylor series, we instead calculate the 
derivative with respect to fi 2 directly. Our calculation confirms the negative sign of 
the leading term on two different volumes and with drastically reduced errors. Second, 
we repeat our fits to Taylor series with vastly increased statistics and simulations at 
additional values of chemical potentials. This allows us to also constrain the 0(fi A ) and 
the sign of the 0(n e ) terms, which we both find to be negative as well. Moreover, our 
leading and next-to-leading order coefficients agree within statistical errors between 
the two methods. Our findings leave little doubt that on coarse N t = 4 lattices with 
staggered fermions the scenario Fig. [1] (right) is realised. 



2 



2 Extracting the fi dependence of the critical point 



Our strategy to determine the /i dependence of the critical quark masses has been 
described in detail in [6]. In order to keep this paper self-contained, we repeat the 
essential formulae. On the lattice, the critical parameters Eqs. (IT|21) get replaced by 
critical couplings, which in turn have a Taylor expansion around the zero density critical 
quark mass, 



(5 c (am, a/i) = (3 c (aml, 0) + cjy (am — amfj) fc (a/i) 2/ , (3) 

k,l=l 

am c (afi) = ami + c ' k ( a ^ 2k ■ ( 4 ) 

k=l 

As in our previous investigations, our observable for criticality is the Binder cumulant 
constructed from the fluctuations of the chiral condensate, 

((5i>il>) 4 ) - - - 

B4 = Mm 2 ' w = ^-w>- (5) 

In the infinite volume limit the Binder cumulant behaves discontinuously, assuming the 
values 1 in a first order regime, 3 in a crossover regime and 1.604 characteristic of the 
3d Ising universality at a second order chiral critical point. In a finite volume the dis- 
continuities are smeared out and flattened, so that B4 passes continuously through the 
critical value, i.e. it is an analytic function of (f3, am, a/i). For a given pair (am, a/i), the 
gauge coupling needs to be fixed to its pseudo-critical value f3 c (am, a/i) by requiring, for 
instance, that the susceptibility ((5X) 2 ) be maximum. We use the equivalent prescrip- 
tion of requiring a vanishing third moment, M3 = ((5X) 3 ) = 0. In the neighbourhood 
of the chiral critical point the Binder cumulant can then be expanded as 

B A (am, a/i) = 1.604 + V" b M (am - am c ) k (a^) 21 , (6) 

k,l=l 

with volume dependent coefficients bki(L). 

For large volumes the approach to the thermodynamic limit is governed by univer- 
sality. Near a critical point the correlation length diverges as £ ~ r~ v , where r is the 
distance to the critical point in the plane of temperature and magnetic field-like vari- 
ables, and v « 0.63 for the 3d Ising universality class. Since (3 is tuned to (3 C always, 
we have r oc \am — am^\ for \x = and r oc |a/i| 2 for am = am^ . B4 is a function of 
the dimensionless ratio L/£, or equivalently (L/£ > ) l ' v . Hence, for the Taylor coefficients 
one expects, for large volumes, the universal scaling behaviour 

b kl (L) = f kl L (k+l ^ , (7) 
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where the fu are independent of the volume. 

It is now straightforward to relate the desired coefficients in Eq. (J4]) to those appear- 
ing in the Binder cumulant Eq. ([6]). The function am c (afi) is defined implicitly by the 
condition that the Binder cumulant stays critical, B 4 (am, a/i) = 1.604. Taking total 
derivatives with respect to /i 2 one obtains 



Ci 



d am c dB A ( dB A \ 1 6 01 

d(afi) 2 d(afi) 2 \dam J b w 



. 1 d 2 am c 1 , . /2 s 

C ' 3 = ^! d\a^) 2 f = ~6^^ 03 + &llC2 + &12C/l + k 1 ^* + 2&2 ° c ' lC2 + &3oC 'i 3 ) • ( 8 ) 

Once these coefficients are known, they need to be converted to the continuum coef- 
ficients Eq. (J2J). This is non-trivial: for fixed temporal lattice extent N t the lattice 
spacings entering the dimensionless am c (n) and am c (0) are different, since the pseudo- 
critical temperature T c (m c (/i), /j) = l/(iV t a(/i)) depends on /j. The relation between 
the continuum and lattice coefficients is thus 

7r 2 c[ 1 dT c (m c (fi), /i) 



Cl 



C 2 



iV 2 amg T c (m§, 0) rf(/i/vrT) 2 ' 
vr 4 c' 2 vr 2 C ; 1 dT c (m c (fi),fi) 
Nj amg ~ ivf amg T c (mg,0) d{n/-nT) 
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1 rf 2 T c K(/i),/i) 
+ 2T C K,0) rf[W7rT) 2 ] 2 • 1 ' 

The pseudo-critical temperature T c (m, /i) can be obtained from /3 c (am,an) by means 
of the two-loop beta-function. Note that for our coarse lattices Nt — 4 this procedure 
is certainly not quantitatively valid, and some non-perturbative beta-function should 
be used. However, we shall see that for the conclusions to be drawn later the two-loop 
beta-function is actually sufficient. 

With these equations, we are ready for numerical evaluations. Our task is to extract 
the coefficients bki of the Binder cumulant expansion, Eq. ([6]), from numerical simula- 
tions, and then convert them to the continuum via Eqs. (I5f§|) . In the following sections 
we present results for two independent ways of measuring the bu- 



3 Calculating derivatives 

When attempting to extract Taylor coefficients by fits of truncated series to data con- 
taining the full /i-dependence at imaginary /i, one has to worry about systematic errors. 
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Fitting to a finite polynomial necessarily distributes the actual /^-dependence across the 
available terms, and the result is only reliable to the extent that the higher terms are 
numerically negligible. In particular, leading order fits run a risk of "averaging" a more 
complicated functional dependence into a single coefficient, especially when the func- 
tion is very flat. These dangers have recently been demonstrated for the pseudo-critical 
coupling in SU(2), where calculations at real /i are possible and can be compared to 
results obtained via analytic continuation from imaginary /t [8]. 

To guard against this, and because the //-dependence of B 4 is extremely weak [6] , we 
now wish to calculate the derivative dB^/d^afi) 2 directly, without recourse to fitting. 
(For dB^/ d(am) this is not necessary, the m-dependence of B 4 is pronounced enough to 
give reliable fits [HE])- In principle, this can be done by expressing the /x— derivative of 
_B 4 analytically through traces of non-local operators, whose expectation values are to 
be evaluated at \i = 0. This technique has been used to calculate the Taylor coefficients 
of the pseudo-critical temperature, the pressure and various susceptibilities in [HI ITU] . 
However, we do not follow this approach here. First, the numerical effort is non-trivial. 
For instance, the trace of the 5 th inverse power of the Dirac operator must be evaluated 
to obtain the derivative of (Sipip) 41 with respect to the quark mass. Moreover, delicate 
cancellations will take place among the rapidly growing number of terms, leading to 
difficulties with optimising the number of noise vectors to be used as accurate stochastic 
estimators of the various traces. 

Instead, we use a new method to estimate the derivative, which turns out to be 
simpler to implement and vastly more efficient. We measure the change AB4 under 
a small variation A(a/t) 2 , thus estimating the ratio of finite differences which, for 
sufficiently small variations, approaches the desired derivative, 



Because the required shift in the couplings is very small, it is adequate and safe to 
use the original Monte Carlo ensemble for amjj, \x = and reweight the results by 
the standard Ferrenberg-Swendsen method [TTJ. Moreover, reweighting to imaginary 
/a the reweighting factors remain real positive and close to 1. With reweighting, the 
fluctuations in the original and the reweighted ensembles are strongly correlated. For 
the calculation of variations this is beneficial. The correlated fluctuations drop out 
of our observable, which now is the change in S 4 , rather than _B 4 itself. As a final 
simplification, we do not calculate the reweighting factors 



AB 4 dB 4 



(10) 



a^ho A(a/i) 2 (9(a/i) 2 




det 



>/4 W/^) 



(11) 
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Figure 2: Finite difference quotients on L = 8 extrapolated to zero as a constant, or 
constant plus (a/i«) 2 correction. The resulting coefficients are listed in Table [U The 
error band corresponds to the highlighted best fits to data at imaginary fi, Table [2j 



exactly, but estimate them through Gaussian noise vectors r] in the standard way, 

p( / i 1 ,^) = <exp(-|r^/ 8 (^)^ /8 (^ih| 2 + l^| 2 )>,- (12) 

As a further check of the method on a more accurately calculable observable, we also 
consider the pseudo-critical coupling /3 C defined by the peak of the plaquette suscep- 
tibility, and calculate <9/3 c /<9(a/i) 2 along the way. Further details and a thorough test 
of this method in the Potts model are described in [12], where the results were found 
to fully agree with the other method of computing derivatives. Since that test the 
technique has also been successfully applied in a similar calculation at finite isospin 
chemical potential [T3] . 

Our calculations for d(3 c / 'd{afi) 2 , dB^j '<9(a/i) 2 were done on 8 3 x 4 and 12 3 x 4 lattices, 
for which we have accumulated about 5 million and 0.5 million trajectories, respectively 
0. The bare quark mass was kept fixed at am = 0.0265, which is close to the critical 
value arriQ and corresponds to m^L ~ 3.4, 5, respectively. (A small offset from am^ 
will affect the //-derivative only via the lowest order mixing term, whose coefficient 
is compatible with zero, as we shall see in the next section.) The results of this 
calculation are shown in Fig. [2j The individual data points correspond to the difference 
quotient reweighted to the given value of /ij. The error bars within one plot are thus 
strongly correlated. The straight lines represent linear extrapolations to zero /ij, and 

1 The 8x4 simulations were performed in about two months on the EGEE computing grid, thanks 
to the support of the CERN IT/GS group. Jobs at different parameter sets were distributed on the 
grid infrastructure where typically ~ 1000 CPUs could be used at the same time for this activity. 
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Coi 


0)2 


rVdof 


&01 


&02 


foi 


/o2 


X 2 /dof 


0.7448(6) 




0.38 


1.36(8) 




0.050(3) 




0.82 


0.7430(6) 


0.31(9) 


0.12 


1.59(12) 


43.4(18.8) 


0.059(4) 


0.059(26) 


0.33 


0.737(3) 




1.10 


2.63(40) 




0.051(8) 




0.50 


0.7430(6) 


0.14(1.10) 


1.37 


3.35(53) 


259.0(148.0) 


0.065(10) 


0.097(56) 


0.44 



Table 1: Taylor coefficients of j3 c (am,a{i) and B^(am, a/i), Eqs. ( I3"l6f7l) . from the ex- 
trapolated finite differences, as shown in Fig. [2j The upper rows are for L = 8, the 
lower rows for L = 12. 

the intercept is the final estimate for the respective derivatives. The extrapolations 
with non-zero slope show the influence of the next-to-leading term. All estimates are 
collected in Table [TJ 

Next, we can check for finite volume scaling between our two volumes. This is done 
in Fig. [3j where the difference quotients are plotted against a volume-rescaled variable. 
The pseudo-critical coupling has a finite thermodynamic limit, so its Taylor coefficients 
should become volume-independent once the lattices are large enough. On the other 
hand, B4 has to scale in a way characteristic of the 3d Ising universality according to 
Eq. ([7]). The rescaled curves in Fig. [3] nearly fall on top of each other, i.e. both large 
volume scaling behaviours are realised to a good approximation between our L = 8 
and L — 12 lattices. 

4 Simulations at imaginary fi 

In this section we revisit our earlier attempt [6] to determine the Taylor coefficients from 
simulations at non- vanishing imaginary /1 = ifii, for which there is no sign problem. In 
this case the data contain no approximation other than the discretisation, and can be 
fitted by the truncated series Eqs. fl3|6|) directly. As in [6], we consider an 8 3 x 4 lattice, 
but we have added some more values of /ij and increased the statistics. We measured 
B± for 42 different pairs (om, o/i) drawn from am £ [0.020, 0.034], /ij £ [0,0.262], each 
pair simulated at 3-5 different /5-values around the pseudo-critical one. Our overall 
statistics is well over 20 million trajectories. 

Fits of Eqs. (13|6|) to various orders are collected in Table [2j The first column specifies 
the fitting range, i.e. the maximal imaginary chemical potential included in the fits. 
We observe that f3 c can be fitted across the entire range, whereas for _B 4 excellent fits 
are obtained up to a/ii = 0.245. Once the fitting range is further increased, the % 2 
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Figure 3: Finite size scaling of the Taylor coefficients between L = 8, 12 lattices. For 
the pseudo-critical coupling, finite size corrections disappear in the thermodynamic 
limit. -B 4 is expected to scale according to Eq. ©. 



abruptly grows, and the values for the coefficients change. Clearly, this indicates the 
need for higher order terms. Furthermore, for large there is a phase transition to an 
unphysical Z(3)-sector [2], which for N t = 4 happens at a/Xj = jj^ ss 0.262. Hence, 
close to this point the scaling behaviour of this Z(3) transition will mix into B4, masking 
that of the transition we are after and leading to large finite volume effects. For these 
reasons we discard fitting ranges a/ij > 0.245. This does not affect f3 c , because it has 
a finite thermodynamic limit across the entire range. 

For the good fits with % 2 /dof < 1.5, all coefficients are stable within errors under 
variations of the fitting range. Interestingly, we have a rather solid signal for a vanishing 
of the lowest order mixing between the m- and /i-dependence in both observables. We 
quote the highlighted fits as our results for the Taylor coefficients from fits to imaginary 
\x data. 

The estimates for the leading terms agree within errors, while the subleading terms 
show some differences. The tendency for both coefficients is towards a smaller absolute 
value than that found by the derivative method. The extended fitting ranges already 
suggest that this is due to the neglect of higher order terms. Nevertheless, within 
one and a half standard deviations, also the sub-leading terms are consistent with 
those from the direct calculation of the derivatives, as is also visible in Fig. [2J We 
conclude that both methods considered here are able to provide the leading coefficient 
of the Taylor series in /1 quantitatively, as well as at least the sign of the sub-leading 
coefficients, though the derivative method is much more economical. 
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| 1 max 


P(arriQ, 0) 


ClO 
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Coi 


C02 


Cll 


X 2 /dof 


0.2 


5.1354(1) 


1.95(3) 


-14.1(2.7) 


-0.75(1) 


0.40(35) 


-0.8(1.2) 


0.53 


0.2 


5.1353(1) 


1.97(2) 


-13.5(2.5) 


-0.75(1) 


0.43(35) 


— 


0.52 


0.2 


5.1352(1) 


1.90(2) 


— 


-0.763(6) 


— 


— 


1.14 


0.245 


5.1354(1) 


1.96(3) 


-14.4(2.1) 


-0.735(7) 


0.84(12) 


-0.31(78) 


0.52 


0.245 


5.1354(1) 


1.97(2) 


-14.5(2.1) 


-0.736(7) 


0.84(11) 


— 


0.51 


0.245 


5.1348(2) 


1.89(2) 


— 


-0.783(5) 


— 


— 


2.17 


0.262 


5.1354(1) 


1.96(3) 


-15.2(2.3) 


-0.729(8) 


1.00(11) 


-0.61(86) 


0.64 


0.262 


5.1354(1) 


1.98(2) 


-15.4(2.3) 


-0.729(7) 


0.99(11) 


— 


0.63 


|/^| max 


am c 


ho 


ho 


boi 




hi 


X 2 /dof 


0.1 


0.0259(4) 


14.0(7) 


_ 


0.70(21) 


_ 


_ 


1.32 


0.2 


0.0257(4) 


14.0(1.5) 


-404.8(119) 


0.38(58) 


-7.2(14.2) 


-63.4(54) 


0.97 


0.2 


0.0259(4) 


14.0(7) 





0.70(21) 








1.32 


0.235 


0.0257(4) 


14.7(1.3) 


-256.5(99) 


0.97(40) 


11.7(6.3) 


-6.9(28) 


1.08 


0.235 


0.0257(4) 


14.9(6) 


-256.4(97) 


0.98(39) 


11.9(6.1) 


— 


1.05 


0.245 


0.0257(4) 


14.6(1.3) 


-255.7(93) 


0.96(38) 


11.4(5.7) 


-8.1(27) 


1.00 


0.245 


0.0257(4) 


14.9(6) 


-254.9(92) 


0.96(37) 


11.6(5.6) 




0.97 


0.255 


0.0255(4) 


15.6(1.5) 


-321.1(118) 


1.46(46) 


22.2(6.8) 


9.4(34) 


1.59 


0.255 


0.0255(4) 


15.3(8) 


-322.4(116) 


1.47(46) 


22.1(6.6) 




1.59 


0.262 


0.0254(5) 


16.5(1.6) 


-311.4(126) 


1.96(47) 


32.2(6.7) 


-33.1(35) 


1.98 


0.262 


0.0252(5) 


15.3(9) 


-308.8(126) 


2.02(48) 


32.4(6.7) 




1.98 



Table 2: Fits of the Taylor expansions /3 c (am,afi), Eq. ([3]), (upper section) and 
B4(am,afi), Eq. ([6]), (lower section) to data at imaginary /i. The first column specifies 
the fitting range. 

5 Comparing and combining approaches 

We now have obtained compatible results from two independent methods of calculat- 
ing the leading Taylor coefficients of power series of critical parameters in afi. We 
stress again that these are not merely two ways of analysing the same data. Rather, 
the simulations required for the two approaches also differ. Furthermore, the deriva- 
tive method Sec. [3] uses small chemical potentials |a/ij| < 0.1, whereas the finite-/^ 
simulations Sec. 4 consider larger values |a/ij| > 0.1. The two approaches are thus 
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ho 
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b 2 


b 3 


hi 


X 2 /dof 


14.8(6) 


-421.7(177) 


49.4(8.3) 


469.6(150) 


-4950.5(4142) 


1.05 


14.9(6) 


-243.8(96) 


47.4(8.2) 


421.9(145) 




1.06 


15.0(7) 


-208.3(184) 


23.5(1.4) 




-1470.7(4495) 


1.33 


15.0(7) 


-156.8(101) 


23.7(1.3) 






1.31 



Table 3: Fits of S 4 (am, a/x), Eq. (jHJ), to imaginary /x data, with the leading coefficients 
arriQ = 0.0259, 6 i = 1-6 fixed to their values determined at /x = 0. 

complementary, and it is reasonable to combine their results in order to maximise the 
available information. To start with, the critical quark mass at zero density, am§, and 
hi are best determined from the /x = calculations of Sec [3l Then, these values can 
be used as input for the fits to a set of imaginary /x data. This is shown in Table [BJ 
where we have fixed omg = 0.0259 and hi = 1.6, according to our best estimates from 
fi = simulations. 

With two coefficients fixed, we can actually consider the next order in _B 4 , i.e. add 
terms C(m 3 , m 2 /x 2 , m/x 4 , /x 6 ), and still get sensible fits. We find 630 = 621 = in all 
combinations of fit parameters. Removing the entirely unconstrained coefficients one 
by one, we observe that 6 3 is the most significant coefficient at this order. Note also 
that allowing for 603, the value of 602 relaxes precisely to the value obtained from the 
derivative method, Table [H The same phenomenon is observed for the /x 2 coefficient 
when releasing /x 4 . Since 601,602,603 ah have the same sign, this is not surprising. At 
imaginary /x = i/Xj consecutive terms alternate in sign, and thus show a tendency to 
cancel each other out. This demonstrates the "averaging" into lower order coefficients 
of higher order terms that get significant but are excluded from a fit. For these reasons 
we believe the highlighted fit in Table El with 603 released, to be the most trustworthy 
one. Consistency between the two methods is also illustrated in Fig. HI where we show 
the data for the reweighted finite difference quotients together with those determined 
from the direct measurements of B4 at finite /x = z/Xj. 

We conclude that by combining both methods discussed here, we are sensitive to the 
/x 6 -term. Its main effect is to give us some confidence in the order of magnitude of the 
/x 4 term determined in our earlier analyses. Moreover, without taking the value of 603 
too seriously, it comes with the same sign as 601 (see Fig. HI), thus leading to a negative 
c' 3 . Hence, the /x 6 -term also contributes to the shrinking of the first order region. 
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Figure 4: Finite difference quotients on L = 8 lattices for am = 0.0265. The left section 
(a/Xj) 2 < 0.01 is determined via reweighting, the right section (ct/ij) 2 > 0.01 via direct 
measurements of at finite \i = ifa. The points at larger values clearly fall below the 
C(/i 4 ) contribution, thus indicating a negative /i 6 -term. 

6 Conclusions 

For the continuum conversion, Eq. (Q, we need the change of temperature at the 
critical point with fi. This is most conveniently extracted by evaluating the pseudo- 
critical couplings at the critical points, fitting as a series in (an) 2 and converting to 
temperatures by the two-loop beta-function. Putting it all together, our final result 
reads§ 

^). 1 _,3,3)(4) 2 -4 7(2 0)(4) 4 -... ( 13 , 

We conclude that our investigation of higher order contributions in fi 2 fully supports 
and sharpens the conclusions of [5J: the light mass range featuring a first order phase 
transition is shrinking with chemical potential. The quoted errors include one standard 
deviation uncertainty on the coefficients of the critical temperature, as well as on the 
parameters mQ,6 i calculated at \i — 0. They do not include, however, the systematic 
error introduced by using the two-loop beta-function. 

Nevertheless, our conclusion is stable against the use of another beta-function, as 
can be seen from Eq. (JHJ). The contribution of the second term to c\ is always negative, 



2 We note that the preliminary result C2 = —12(6) reported earlier at various conferences [15j 
corresponds to the value obtained from fits to imaginary \x data alone. 
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irrespective of the beta-function used. For C2, the first term is numerically dominant 
and the second is guaranteed to be negative too, no matter what the beta-function. 
The third term presently also gives a negative contribution, and it is hard to see how 
a change in beta-function would change the sign and overhaul the contribution of the 
first two terms. Finally, we have good evidence that c' 3 is also negative, i.e. the /x 6 
term further shrinks the first order region, at least in lattice unitsjj Thus, all terms 
appear to have negative signs, confirming our conclusion irrespective of their precise 
quantitative values. 

Our result therefore demonstrates beyond reasonable doubt that for Nf = 3 the 
region of first order chiral phase transitions shrinks for small to moderate fi/T, as in 
Fig. [H (right), and hence there is no chiral critical point at chemical potentials [i<T 
in a theory with three flavours of staggered fermions and m > m% on an N t = 4 lattice. 
In order to draw conclusions for physical QCD, two further steps are necessary. Firstly, 
the investigation has to be generalised to Nf — 2 + 1. Secondly, and most importantly, 
it has to be repeated on finer lattices, since N t = 4 with a ~ 0.3 fm is too coarse to 
accurately reflect continuum physics. It is by no means excluded that cut-off effects 
are stronger than finite density effects on coarse lattices, and hence a change of sign in 
the curvature of the chiral critical surface in the continuum is a possibility. Simulations 
to address these issues are in progress. Finally, we emphasise that our findings only 
concern a chiral critical point, i.e. a point on the critical surface in Fig. [T] bounding 
the region of first-order chiral transitions. Our results do not exclude an additional 
critical surface in Fig. [T] (right), not analytically connected to the region of chiral phase 
transitions, which could cause a critical point and phase transition of a different nature 
to exist. 
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